!Résolution de l'équation intégrale linéaire de solution \exp(-x). Pas de recollement.
!Equation utilisée: de type "grands c"
module fonctions
implicit none
integer, parameter :: points = 3000 !nombre de points d'évaluation des fonctions
integer, parameter :: maxit = 80 !nombre maximal d'itérations
double precision, parameter :: infty = 20.d0 !valeur numérique pour l'infini
double precision, parameter :: lambda = 1.00d0 !vitesse de convergence (=1.0) versus stabilité (=0.01...) (?)
contains

!------------------------------------------------------------------------
!Méthode de simpson pour calculer une intégrale connaissant uniquement
!quelques points discrets. Explication des variables:
!  y       : tableau de dimension (n+1) contenant les valeurs y(xi)
!  h       : largeur du pas d'intégration
!  n       : nombre de pas d'intégration
!  somme   : valeur de l'intégrale sur la base des points fournis
subroutine simptab(y,h,n,somme)
implicit none
double precision, intent(in), dimension(0:points) :: y
double precision, intent(in) :: h
integer, intent(in) :: n
double precision, intent(out) :: somme
integer :: i
double precision :: wgt
somme = y(0)
somme = somme + y(n)
if (n==0) somme = 0.d0
if (n==1) somme = h*somme/2.d0
wgt = 2.d0
if (n>1) then
	do i=1,n-1
	  if (wgt == 2.d0) then
	    wgt = 4.d0
	  else
	    wgt = 2.d0
	  end if
	  somme = somme + wgt*y(i)
	end do
	somme = h*somme/3.d0
end if
return
end subroutine simptab

!------------------------------------------------------------------------
!Méthode du trapèze pour calculer une intégrale connaissant uniquement
!quelques points discrets. Explication des variables:
!  y       : tableau de dimension (n+1) contenant les valeurs y(xi)
!  h       : largeur du pas d'intégration
!  n       : nombre de pas d'intégration
!  somme   : valeur de l'intégrale sur la base des points fournis
subroutine traptab(y,h,n,somme)
implicit none
double precision, intent(in), dimension(0:points) :: y
double precision, intent(in) :: h
integer, intent(in) :: n
double precision, intent(out) :: somme
integer :: i
somme = y(0)
somme = somme + y(n)
if (n==0) somme = 0.d0
if (n > 1) then
	do i=1,n-1
	  somme = somme + 2.d0*y(i)
	end do
end if
somme = h*somme/2.d0
return
end subroutine traptab

end module fonctions

!----------------------------------------------------------------------------------
!----------------------------------------------------------------------------------
!résolution de f(x) = 1- \int_0^x dy f(y), solution: f(x) = \exp(-x)

program anneqinttest
use fonctions
implicit none
double precision, dimension(0:points) :: fout,fin
double precision :: h,err,err2
integer :: i,j

h = infty/dble(points)
fin(i) = 1.d0 !condition initiale
do j=1,maxit
	do i=0,points
		call traptab(fin,h,i,fout(i))
	end do
	fout = 1.d0 - fout
	err = sum(dabs(fin-fout))/dabs(sum(fout))
	fin = (1.d0 - lambda)*fin + lambda*fout
	err2 = 0.d0 	!erreur à la solution exponentielle
	do i=0,points
		err2 = err2 + dabs(fin(i) - dexp(-h*dble(i)))
	end do
	open(unit=1, file='ann-eqinttest.txt')
	do i=0,points
		write(1,*) h*dble(i),fout(i)
	end do
	close(1)
	write(*,'(" Iteration: ",i5,", err.relative: ",e9.2,", |f_n(x) - exp(-x)|/points: ",e9.2)') j,err,err2/dble(points)
end do

write(*,*) "Fin du programme."
read(*,*)
end program anneqinttest